﻿using System;
using System.Collections.Generic;
using LatoolNet;

namespace NumericalSolution {
  class Poisson1DFEM {

    private int m_numberOfElements;
    private int m_numberOfNodes;
    private int m_numberOfNodesInAnElement;
    private NumericalIntegration m_integration;
    private double m_beta = Double.NaN;
    private Element[] m_elements;
    private Dictionary<int, Node> m_nodeMap = new Dictionary<int, Node>();
    private Matrix m_A;
    private Matrix m_b;

    public int NumberOfElements {
      get { return m_numberOfElements; }
      set {
        m_numberOfElements = value;
        m_numberOfNodes = value + 1;
      }
    }

    public int NumberOfNodes {
      get { return m_numberOfNodes; }
    }

    public int NumberOfNodesInAnElement {
      get { return m_numberOfNodesInAnElement; }
      set { m_numberOfNodesInAnElement = value; }
    }

    public double Beta {
      get { return m_beta; }
      set { m_beta = value; }
    }

    public NumericalIntegration NumericalIntegration {
      get { return m_integration; }
      set { m_integration = value; }

    }
    public void CreateMesh() {

      double x_begin = 0.0;
      double x_end = 1.0;
      double x_delta = (x_end - x_begin) / m_numberOfElements;

      m_elements = new Element[m_numberOfElements];
      int number = 0;
      for (int i = 0; i < m_numberOfElements; i++) {
        Element e = new Element(m_numberOfNodesInAnElement);
        m_elements[i] = e;
        e.Number = number++;
        e.Node[0].X = i * x_delta;
        e.Node[1].X = (i + 1) * x_delta;

        if (m_nodeMap.ContainsKey(q(e.Number,0)) == false) {
          m_nodeMap.Add(q(e.Number, 0), e.Node[0]);
        }
        if (m_nodeMap.ContainsKey(q(e.Number,1)) == false) {
          m_nodeMap.Add(q(e.Number, 1), e.Node[1]);
        }
      }
    }

    public void AssembleMatrix() {

      m_A = new Matrix(m_numberOfNodes, m_numberOfNodes, 3);
      m_b = new Matrix(m_numberOfNodes, 1);

      for (int i = 0; i < m_numberOfElements; i++) {
        Matrix A_e = new Matrix(m_numberOfNodesInAnElement, 
                                m_numberOfNodesInAnElement);
        Matrix b_e = new Matrix(m_numberOfNodesInAnElement, 1);

        Element element = m_elements[i];
        AssembleLocalMatrix(A_e, b_e, element);
        AssembleGlobalMatrix(A_e, b_e, element);
      }

    }

    private void AssembleLocalMatrix(Matrix A_e, 
                                     Matrix b_e, 
                                     Element element) {

      for (int integrationIndex = 0;
           integrationIndex < m_integration.NumberOfIntegrationPoints;
           integrationIndex++) {

        Integrands(A_e, b_e, element,
                   m_integration.IntegrationPoint(integrationIndex),
                   m_integration.Weight(integrationIndex));
      }
    }

    private void Integrands(Matrix A_e, Matrix b_e, Element element,
                            double xi, double weight) {
      double h = element.Node[1].X - element.Node[0].X;
      double J = h / 2;

      for (int r = 0; r < m_numberOfNodesInAnElement; r++) {
        for (int s = 0; s < m_numberOfNodesInAnElement; s++) {
          A_e[r, s] += dPhi(r, xi) * (2 / h) *
                       dPhi(s, xi) * (2 / h) *
                       J * weight;
        }
        b_e[r, 0] += m_beta * Phi(r, xi) * J * weight;
      }
    }

    private void AssembleGlobalMatrix(Matrix A_e, Matrix b_e, Element e) {

      if (e.Number == 0) {
        A_e[0, 0] = 1;
        A_e[0, 1] = 0;
        b_e[0, 0] = 0;
      }
      if (e.Number == m_numberOfElements - 1) {
        A_e[1, 1] = 1;
        A_e[1, 0] = 0;
        b_e[1, 0] = 1;
      }
      for (int r = 0; r < m_numberOfNodesInAnElement; r++) {
        for (int s = 0; s < m_numberOfNodesInAnElement; s++) {
          m_A[q(e.Number, r), q(e.Number, s)] += A_e[r, s];
        }
        m_b[q(e.Number, r), 0] += b_e[r, 0];
      }
    }

    public Node[] Solve() {

      //LUFactorization.Solve(m_A, m_b);
      LinearEquation.Solve(m_A, m_b);

      List<Node> nodelist = new List<Node>();
      for (int i = 0; i < m_numberOfNodes; i++) {
        Node node = m_nodeMap[i];
        node.Value = m_b[i, 0];
        nodelist.Add(node);
      }

      m_A.Dispose();
      m_b.Dispose();
      return nodelist.ToArray();
    }

    private double Phi(int type, double xi) {
      if (type == 0) {
        return 0.5 * (1 - xi);
      } else if (type == 1) {
        return 0.5 * (1 + xi);
      }
      throw new ArgumentOutOfRangeException();
    }

    private double dPhi(int type, double xi) {
      if (type == 0) {
        return -0.5;
      } else if (type == 1) {
        return 0.5;
      }
      throw new ArgumentOutOfRangeException();
    }

    private int q(int elementNumber, int localNumber) {
      return elementNumber + localNumber;
    }
  }
}
